# %%
from typing import List
import dill
import numpy as np
import scipy as sc
from EF2D import *
from EF2D.tools import *
from icecream import ic
# l:/FEACODE/DevEnv/PortablePython310/App/Python/python.exe -m tests.ElemTests.c3d8.my_c3d8
# 其他数据文件在:
initNumpyPrintOptions()
nodes, elements, nsets, elsets,surfaces=read_abaqus_inp("L:/EF2D/tests/ElemTests/c3d8/250407/getElemMatrix - 副本.inp")
node_dict={}
for n in nodes:
    node_dict[n[0]]=Node(label=n[0],
                        x=n[1],y=n[2],z=n[3])
globaldof=nodes.shape[0]*3

elem_list:List[C3D8]=[]
for e in elements:
    E=C3D8(label=e[0],
         node1=node_dict[e[1]],
         node2=node_dict[e[2]],
         node3=node_dict[e[3]],
         node4=node_dict[e[4]],
         node5=node_dict[e[5]],
         node6=node_dict[e[6]],
         node7=node_dict[e[7]],
         node8=node_dict[e[8]],
         E=77200,nu=0.3)
    elem_list.append(E)
ic(len(elem_list))
elem=elem_list[0]

# 装配刚度矩阵
Kg=np.zeros((globaldof,globaldof))
Kg[np.ix_(elem.ElemDofIndexs,elem.ElemDofIndexs)]+=elem.Ke_bbar

# 节点载荷10N
load_dict={get_ithdof_index(i,3,2)+1:10.0 for i in [2,4,6,8]}
# fix 四个节点的x,y,z方向的位移
bc_dict={ get_ithdof_index(i,3,j)+1:0.0 for i in [1,3,5,7] for j in [1,2,3]}

u,f=Method2(Kg,bc_dict,load_dict)
U1=u[::3]
U2=u[1::3]
U3=u[2::3]
# ic(U1,U2,U3)
plot1d(U2)
# 
# imShowMat(elem.Ke,"no b-bar ke")
sc.io.savemat("tests/ElemTests/c3d8/250407/my_c3d8_element_stiff250407.mat",
              {'Ke':elem.Ke,"Ke_bbar":elem.Ke_bbar,'my_u':u,'my_f':f,
               "globalKg":Kg,"elemIndex":elem.ElemDofIndexs})

# %%

u_,f_=Method2(elem.Ke_bbar,bc_dict,load_dict)

U1_=u_[::3]
U2_=u_[1::3]
U3_=u_[2::3]

plot1d(U2_)
# sc.io.savemat("tests/ElemTests/c3d8/bbar_my_c3d8_element_stiffness.mat",{'Ke_bbar':elem.Ke_bbar})

# %%
pp=sc.io.loadmat("L:/EF2D/tests/ElemTests/c3d8/abaqus_ElementStiffness.mat")
abaqus_ke=pp["Element1Stiffness"]
imShowMat(elem.Ke_bbar-abaqus_ke,"different of my-bbar c3d8 and abaqus c3d8")
# %%
imShowMat(elem.Ke-abaqus_ke,"different of my-bbar c3d8 and abaqus c3d8")


# %%
imShowMat(np.abs(elem.Ke-elem.Ke_bbar),"elem.Ke-elem.Ke_bbar")
